#include"erpa.h"
#include"constants.h"
#include"help.h"
ERPA::ERPA(const System_Tz_HF* Ptr):pSystem_Tz_HF(Ptr)
{
  vector<TwoBodyState_ph> temp;
  Jmax=0;
  while(1)
    {
      if(set_ph_same(Jmax,0,temp)>0||set_ph_same(Jmax,1,temp)>0
	 ||set_ph_np(Jmax,0,temp)>0||set_ph_np(Jmax,1,temp)>0
	 ||set_ph_pn(Jmax,0,temp)>0||set_ph_pn(Jmax,1,temp)>0
	 )
	++Jmax;
      else
	break;
    }
  ph_sames.resize(2*Jmax);
  ph_nps.resize(2*Jmax);
  ph_pns.resize(2*Jmax);

  A_sames.resize(2*Jmax);
  B_sames.resize(2*Jmax);
  X_sames.resize(2*Jmax);
  Y_sames.resize(2*Jmax);
  Omega_sames.resize(2*Jmax);

  Anpnps.resize(2*Jmax);
  Apnpns.resize(2*Jmax);
  Bnppns.resize(2*Jmax);
  T_plus_Xs.resize(2*Jmax);
  T_plus_Ys.resize(2*Jmax);
  T_plus_Omegas.resize(2*Jmax);
  T_minus_Xs.resize(2*Jmax);
  T_minus_Ys.resize(2*Jmax);
  T_minus_Omegas.resize(2*Jmax);

  G_sames.resize(2*Jmax);
  G_nps.resize(2*Jmax);
  G_pns.resize(2*Jmax);
}

int ERPA::set_ph_same(int J,int Par,vector<TwoBodyState_ph>&ph_same)
{
  const  HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;
  ph_same.clear();
  for(int h=0;h<=pSystem_Tz_HF->FermiSurface;++h)
    {
      for(int p=pSystem_Tz_HF->FermiSurface+1;p<HFOrbs.size();++p)
	{
      	  if(HFOrbs[h].tz!=HFOrbs[p].tz) continue;
	  if((HFOrbs[h].l+HFOrbs[p].l-Par)%2!=0) continue;
	  if(abs(HFOrbs[h].jj-HFOrbs[p].jj)>2*J) continue;
	  if((HFOrbs[h].jj+HFOrbs[p].jj)<2*J) continue;
	  TwoBodyState_ph temp(p,h);
	  ph_same.push_back(temp);
	}
    }
  return ph_same.size();
}
int ERPA::set_ph_np(int J,int Par,vector<TwoBodyState_ph>&ph_np)
{
  const HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;
  ph_np.clear();
  for(int h=0;h<=pSystem_Tz_HF->FermiSurface;++h)
    {
      for(int p=pSystem_Tz_HF->FermiSurface+1;p<HFOrbs.size();++p)
	{
	  //tz=-1 for proton,1 for neutron
	  if((HFOrbs[h].tz!=-1)||(HFOrbs[p].tz!=1)) continue;
	  if((HFOrbs[h].l+HFOrbs[p].l-Par)%2!=0) continue;
	  if(abs(HFOrbs[h].jj-HFOrbs[p].jj)>2*J) continue;
	  if((HFOrbs[h].jj+HFOrbs[p].jj)<2*J) continue;
	  TwoBodyState_ph temp(p,h);
	  ph_np.push_back(temp);
	}
    }
  return ph_np.size();
}
int ERPA::set_ph_pn(int J,int Par,vector<TwoBodyState_ph>&ph_pn)
{
  const HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;
  ph_pn.clear();
  for(int h=0;h<=pSystem_Tz_HF->FermiSurface;++h)
    {
      for(int p=pSystem_Tz_HF->FermiSurface+1;p<HFOrbs.size();++p)
	{
	  //tz=-1 for proton,1 for neutron
	  if((HFOrbs[h].tz!=1)||(HFOrbs[p].tz!=-1)) continue;
	  if((HFOrbs[h].l+HFOrbs[p].l-Par)%2!=0) continue;
	  if(abs(HFOrbs[h].jj-HFOrbs[p].jj)>2*J) continue;
	  if((HFOrbs[h].jj+HFOrbs[p].jj)<2*J) continue;
	  TwoBodyState_ph temp(p,h);
	  ph_pn.push_back(temp);
	}
    }
  return ph_pn.size();
}

double ERPA::e(int bra,int ket)
{
  const HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;
  if( !HFOrbs[bra].isInAGroupWith(HFOrbs[ket]) )
    return 0;
    
  double t_temp=pSystem_Tz_HF->get1B(bra,ket);
  double v_temp=0;
  for(int i=0;i<HFOrbs.size();i++)
    {
      for(int k=0;k<HFOrbs.size();k++)
	{
	  if( !HFOrbs[i].isInAGroupWith(HFOrbs[k]) ) continue;
	  int j_min=abs(HFOrbs[ket].jj-HFOrbs[i].jj)/2;
	  int j_max=(HFOrbs[ket].jj+HFOrbs[i].jj)/2;
	  for(int j=j_min;j<=j_max;j++)
	    v_temp+=pSystem_Tz_HF->get2B(bra,i,ket,k,j)*density_matrix(i,k)*(2*j+1);
	}
    }
  v_temp/=(HFOrbs[bra].jj+1);
  return t_temp+v_temp;
}
void ERPA::set_e()
{
  const HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;
#pragma omp parallel
  {
#pragma omp for schedule(dynamic)
    for(int i=0;i<=pSystem_Tz_HF->FermiSurface;i++)
      {
	for(int j=i;j<=pSystem_Tz_HF->FermiSurface;j++)
	  {
	    e_h(i,j)=e(i,j);
	    e_h(j,i)=e_h(i,j);
	  }
      }
#pragma omp for schedule(dynamic)
    for(int i=pSystem_Tz_HF->FermiSurface+1;i<HFOrbs.size();i++)
      {
	for(int j=i;j<HFOrbs.size();j++)
	  {
	    e_p(i-pSystem_Tz_HF->FermiSurface-1,j-pSystem_Tz_HF->FermiSurface-1)=e(i,j);
	    e_p(j-pSystem_Tz_HF->FermiSurface-1,i-pSystem_Tz_HF->FermiSurface-1)=e_p(i-pSystem_Tz_HF->FermiSurface-1,j-pSystem_Tz_HF->FermiSurface-1);
	  }
      }
  }
}
double ERPA::a(int bra_p,int bra_h,int ket_p,int ket_h,int J)
{
  const HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;
  double temp=density_matrix(bra_h,ket_h)*e_p(bra_p-pSystem_Tz_HF->FermiSurface-1,ket_p-pSystem_Tz_HF->FermiSurface-1)+density_matrix(bra_p,ket_p)*e_h(bra_h,ket_h);

  double tempd=0;
  if(bra_h==ket_h&&(HFOrbs[bra_p].isInAGroupWith(HFOrbs[ket_p])))
    {
      for(int i=0;i<p_ljjtzs[table[bra_p].block_number].ord_list.size();i++)
	{
	  tempd+=e_p(p_ljjtzs[table[bra_p].block_number].ord_list[i]-pSystem_Tz_HF->FermiSurface-1,ket_p-pSystem_Tz_HF->FermiSurface-1)*density_matrix(bra_p,p_ljjtzs[table[bra_p].block_number].ord_list[i]);
	  tempd+=e_p(p_ljjtzs[table[bra_p].block_number].ord_list[i]-pSystem_Tz_HF->FermiSurface-1,bra_p-pSystem_Tz_HF->FermiSurface-1)*density_matrix(ket_p,p_ljjtzs[table[bra_p].block_number].ord_list[i]);
	}
    }
  temp-=0.5*tempd;
  tempd=0;
  if(bra_p==ket_p&&(table[bra_h].block_number==table[ket_h].block_number))
    {
      for(int i=0;i<h_ljjtzs[table[bra_h].block_number].ord_list.size();i++)
	{
	  tempd+=e_h(h_ljjtzs[table[bra_h].block_number].ord_list[i],ket_h)*density_matrix(bra_h,h_ljjtzs[table[bra_h].block_number].ord_list[i]);
	  tempd+=e_h(h_ljjtzs[table[bra_h].block_number].ord_list[i],bra_h)*density_matrix(ket_h,h_ljjtzs[table[bra_h].block_number].ord_list[i]);
	}
    }
  temp-=0.5*tempd;
  
  int bra_p_jj=HFOrbs[bra_p].jj;int bra_h_jj=HFOrbs[bra_h].jj;
  int ket_p_jj=HFOrbs[ket_p].jj;int ket_h_jj=HFOrbs[ket_h].jj;
  int j_min=max(abs(ket_p_jj-bra_h_jj),abs(bra_p_jj-ket_h_jj))/2;
  int j_max=min(ket_p_jj+bra_h_jj,bra_p_jj+ket_h_jj)/2;
  for(int k=j_min;k<=j_max;++k)
    {
      int phase=((ket_p_jj+bra_h_jj)/2-k)%2?-1:1;
      for(int i=0;i<h_ljjtzs[table[bra_h].block_number].ord_list.size();i++)
	for(int j=0;j<h_ljjtzs[table[ket_h].block_number].ord_list.size();j++)
	  {
	    temp+=phase*(2*k+1)*gsl_sf_coupling_6j(bra_p_jj,bra_h_jj,2*J,ket_p_jj,ket_h_jj,2*k)*pSystem_Tz_HF->get2B(bra_p,h_ljjtzs[table[ket_h].block_number].ord_list[j],h_ljjtzs[table[bra_h].block_number].ord_list[i],ket_p,k)*density_matrix(h_ljjtzs[table[ket_h].block_number].ord_list[j],ket_h)*density_matrix(h_ljjtzs[table[bra_h].block_number].ord_list[i],bra_h);
	  }
      
      for(int i=0;i<p_ljjtzs[table[ket_p].block_number].ord_list.size();i++)
	for(int j=0;j<p_ljjtzs[table[bra_p].block_number].ord_list.size();j++)
	  {
	    temp+=phase*(2*k+1)*gsl_sf_coupling_6j(bra_p_jj,bra_h_jj,2*J,ket_p_jj,ket_h_jj,2*k)*pSystem_Tz_HF->get2B(p_ljjtzs[table[bra_p].block_number].ord_list[j],ket_h,bra_h,p_ljjtzs[table[ket_p].block_number].ord_list[i],k)*density_matrix(p_ljjtzs[table[ket_p].block_number].ord_list[i],ket_p)*density_matrix(p_ljjtzs[table[bra_p].block_number].ord_list[j],bra_p);
	  }
      
      for(int i=0;i<p_ljjtzs[table[bra_p].block_number].ord_list.size();i++)
	for(int j=0;j<h_ljjtzs[table[ket_h].block_number].ord_list.size();j++)
	  {
	    temp-=phase*(2*k+1)*gsl_sf_coupling_6j(bra_p_jj,bra_h_jj,2*J,ket_p_jj,ket_h_jj,2*k)*pSystem_Tz_HF->get2B(p_ljjtzs[table[bra_p].block_number].ord_list[i],h_ljjtzs[table[ket_h].block_number].ord_list[j],bra_h,ket_p,k)*density_matrix(bra_p,p_ljjtzs[table[bra_p].block_number].ord_list[i])*density_matrix(ket_h,h_ljjtzs[table[ket_h].block_number].ord_list[j]);
	  }
      
      for(int i=0;i<p_ljjtzs[table[ket_p].block_number].ord_list.size();i++)
	for(int j=0;j<h_ljjtzs[table[bra_h].block_number].ord_list.size();j++)
	  {
	    temp-=phase*(2*k+1)*gsl_sf_coupling_6j(bra_p_jj,bra_h_jj,2*J,ket_p_jj,ket_h_jj,2*k)*pSystem_Tz_HF->get2B(bra_p,ket_h,h_ljjtzs[table[bra_h].block_number].ord_list[j],p_ljjtzs[table[ket_p].block_number].ord_list[i],k)*density_matrix(ket_p,p_ljjtzs[table[ket_p].block_number].ord_list[i])*density_matrix(bra_h,h_ljjtzs[table[bra_h].block_number].ord_list[j]);
	  }
    }
  return temp;
}

double ERPA::b(int bra_p,int bra_h,int ket_p,int ket_h,int J)
{
  const HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;
  double temp=0.0;
  int bra_p_jj=HFOrbs[bra_p].jj;int bra_h_jj=HFOrbs[bra_h].jj;
  int ket_p_jj=HFOrbs[ket_p].jj;int ket_h_jj=HFOrbs[ket_h].jj;
  int j_min=max(abs(ket_p_jj-bra_p_jj),abs(bra_h_jj-ket_h_jj))/2;
  int j_max=min(ket_p_jj+bra_p_jj,bra_h_jj+ket_h_jj)/2;
  for(int k=j_min;k<=j_max;++k)
    {	 
      int phase=((ket_p_jj+bra_h_jj)/2+k+J)%2?-1:1;
      for(int i=0;i<h_ljjtzs[table[ket_h].block_number].ord_list.size();i++)
	for(int j=0;j<h_ljjtzs[table[bra_h].block_number].ord_list.size();j++)
	  {
	    temp+=phase*(2*k+1)*gsl_sf_coupling_6j(bra_p_jj,bra_h_jj,2*J,ket_h_jj,ket_p_jj,2*k)*pSystem_Tz_HF->get2B(bra_p,ket_p,h_ljjtzs[table[bra_h].block_number].ord_list[j],h_ljjtzs[table[ket_h].block_number].ord_list[i],k)*density_matrix(h_ljjtzs[table[ket_h].block_number].ord_list[i],ket_h)*density_matrix(h_ljjtzs[table[bra_h].block_number].ord_list[j],bra_h);
	  }
      
      for(int i=0;i<p_ljjtzs[table[ket_p].block_number].ord_list.size();i++)
	for(int j=0;j<p_ljjtzs[table[bra_p].block_number].ord_list.size();j++)
	  {
	    temp+=phase*(2*k+1)*gsl_sf_coupling_6j(bra_p_jj,bra_h_jj,2*J,ket_h_jj,ket_p_jj,2*k)*pSystem_Tz_HF->get2B(p_ljjtzs[table[bra_p].block_number].ord_list[j],p_ljjtzs[table[ket_p].block_number].ord_list[i],bra_h,ket_h,k)*density_matrix(p_ljjtzs[table[ket_p].block_number].ord_list[i],ket_p)*density_matrix(p_ljjtzs[table[bra_p].block_number].ord_list[j],bra_p);
	  }
      
      for(int i=0;i<h_ljjtzs[table[ket_h].block_number].ord_list.size();i++)
	for(int j=0;j<p_ljjtzs[table[bra_p].block_number].ord_list.size();j++)
	  {
	    temp-=phase*(2*k+1)*gsl_sf_coupling_6j(bra_p_jj,bra_h_jj,2*J,ket_h_jj,ket_p_jj,2*k)*pSystem_Tz_HF->get2B(p_ljjtzs[table[bra_p].block_number].ord_list[j],ket_p,bra_h,h_ljjtzs[table[ket_h].block_number].ord_list[i],k)*density_matrix(h_ljjtzs[table[ket_h].block_number].ord_list[i],ket_h)*density_matrix(p_ljjtzs[table[bra_p].block_number].ord_list[j],bra_p);
	  }

      for(int i=0;i<h_ljjtzs[table[bra_h].block_number].ord_list.size();i++)
	for(int j=0;j<p_ljjtzs[table[ket_p].block_number].ord_list.size();j++)
	  {
	    temp-=phase*(2*k+1)*gsl_sf_coupling_6j(bra_p_jj,bra_h_jj,2*J,ket_h_jj,ket_p_jj,2*k)*pSystem_Tz_HF->get2B(bra_p,p_ljjtzs[table[ket_p].block_number].ord_list[j],h_ljjtzs[table[bra_h].block_number].ord_list[i],ket_h,k)*density_matrix(h_ljjtzs[table[bra_h].block_number].ord_list[i],bra_h)*density_matrix(p_ljjtzs[table[ket_p].block_number].ord_list[j],ket_p);
	  }
    }
  return temp;
}

double ERPA::density_matrix(int bra,int ket)
{
  const HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;
  if(bra<=pSystem_Tz_HF->FermiSurface)
    {
      if(ket>pSystem_Tz_HF->FermiSurface)
	return 0;
      else
	{
	  if(table[bra].block_number!=table[ket].block_number)
	    return 0;
	  else
	    return density_matrix_h[table[bra].block_number](table[bra].n,table[ket].n);
	}
    }
  if(bra>pSystem_Tz_HF->FermiSurface)
    {
      if(ket<=pSystem_Tz_HF->FermiSurface)
	{
	  return 0;
	}
      else
	{
	  if(table[bra].block_number!=table[ket].block_number)
	    return 0;
	  else
	    return density_matrix_p[table[bra].block_number](table[bra].n,table[ket].n);
	}
    }
}
void ERPA::set_AB_same(int J,int par)
{
  set_ph_same(J,par,ph_sames[2*J+par]);
  int dim=ph_sames[2*J+par].size();
  A_sames[2*J+par].resize(dim,dim);
  B_sames[2*J+par].resize(dim,dim);
  G_sames[2*J+par].resize(dim,dim);

  for(int i=0;i<dim;i++)
    {
      for(int j=i;j<dim;j++)
	{
	  int bra_p=ph_sames[2*J+par][i].a;
	  int bra_h=ph_sames[2*J+par][i].b;
	  int ket_p=ph_sames[2*J+par][j].a;
	  int ket_h=ph_sames[2*J+par][j].b;
	  G_sames[2*J+par](i,j)=(((bra_p==ket_p)?density_matrix(bra_h,ket_h):0)-((bra_h==ket_h)?density_matrix(bra_p,ket_p):0));
	  G_sames[2*J+par](j,i)=G_sames[2*J+par](i,j);
	}
    }
  
#pragma omp parallel for schedule(dynamic) collapse(2) 
  for(int A_OR_B=0;A_OR_B<2;A_OR_B++)
    {
      for(int i=0;i<dim;++i)
	{
	  for(int j=i;j<dim;++j)
	    {
	      int bra_p=ph_sames[2*J+par][i].a;
	      int bra_h=ph_sames[2*J+par][i].b;
	      int ket_p=ph_sames[2*J+par][j].a;
	      int ket_h=ph_sames[2*J+par][j].b;
	      if(A_OR_B)
		{
		  A_sames[2*J+par](i,j)=a(bra_p,bra_h,ket_p,ket_h,J);
		  A_sames[2*J+par](j,i)=A_sames[2*J+par](i,j);
		}
	      else
		{
		  B_sames[2*J+par](i,j)=b(bra_p,bra_h,ket_p,ket_h,J);
		  B_sames[2*J+par](j,i)=B_sames[2*J+par](i,j);
		}
	    }
	}
    }
}

bool ERPA::cal_same(int J,int par)
{
  set_AB_same(J,par);
  if(ph_sames[2*J+par].empty())
    {
      A_sames[2*J+par].resize(0,0);
      B_sames[2*J+par].resize(0,0);
      X_sames[2*J+par].resize(0,0);
      Y_sames[2*J+par].resize(0,0);
      Omega_sames[2*J+par].resize(0,1);
      return 0;
    }

  SelfAdjointEigenSolver<MatrixXd> SA_ES(G_sames[2*J+par]);
  MatrixXd S=SA_ES.eigenvectors();
  VectorXd Lambda=SA_ES.eigenvalues();

  A_sames[2*J+par]=S.transpose()*A_sames[2*J+par]*S;
  B_sames[2*J+par]=S.transpose()*B_sames[2*J+par]*S;
  for(int i=0;i<Lambda.size();i++)
    {
      A_sames[2*J+par].row(i)*=1.0/sqrt(Lambda(i));
      A_sames[2*J+par].col(i)*=1.0/sqrt(Lambda(i));
      B_sames[2*J+par].row(i)*=1.0/sqrt(Lambda(i));
      B_sames[2*J+par].col(i)*=1.0/sqrt(Lambda(i));
    }
  
  int dim=ph_sames[2*J+par].size();
  MatrixXcd R(2*dim,2*dim);
  for(int i=0;i<dim;++i)
    {
      for(int j=0;j<dim;++j)
	{
	  R(i,j)=complex<double>(A_sames[2*J+par](i,j),0);
	  R(dim+i,dim+j)=complex<double>(-A_sames[2*J+par](i,j),0);
	  R(i,dim+j)=complex<double>(B_sames[2*J+par](i,j),0);
	  R(dim+i,j)=complex<double>(-B_sames[2*J+par](i,j),0);
	}
    }
  ComplexEigenSolver<MatrixXcd> es(R);
  VectorXcd eigen_values=es.eigenvalues();
  MatrixXcd eigen_vectors=es.eigenvectors();
  Omega_sames[2*J+par].resize(dim,1);
  X_sames[2*J+par].resize(dim,dim);
  Y_sames[2*J+par].resize(dim,dim);
  int j=0;int zero_count=0;
  for(int i=0;i<2*dim;++i)
    {
      double norm=eigen_vectors.col(i).head(dim).real().squaredNorm()-eigen_vectors.col(i).tail(dim).real().squaredNorm();
      if(norm>0)
      	{
      	  X_sames[2*J+par].col(j)=eigen_vectors.col(i).head(dim).real();
      	  Y_sames[2*J+par].col(j)=eigen_vectors.col(i).tail(dim).real();
      	  X_sames[2*J+par].col(j)/=sqrt(norm);
      	  Y_sames[2*J+par].col(j)/=sqrt(norm);
      	  Omega_sames[2*J+par](j)=eigen_values(i).real();
      	  j++;
      	}


      //to handle the spurious state of negative 1 caused by the breaking symmetry of translation.
      // if(abs(eigen_values(i).real())<3)
      // 	{
      // 	  if(zero_count)
      // 	    continue;
      // 	  else
      // 	    {
      // 	      ++zero_count;
      // 	      Omega_sames[2*J+par](j)=eigen_values(i).real();
      // 	      X_sames[2*J+par].col(j)=eigen_vectors.col(i).head(dim).real();
      // 	      Y_sames[2*J+par].col(j)=eigen_vectors.col(i).tail(dim).real();
      // 	      double norm=X_sames[2*J+par].col(j).squaredNorm()-Y_sames[2*J+par].col(j).squaredNorm();
      // 	      if(norm>0)
      // 		{
      // 		  X_sames[2*J+par].col(j)/=sqrt(norm);
      // 		  Y_sames[2*J+par].col(j)/=sqrt(norm);
      // 		}
      // 	      else
      // 		{
      // 		  Omega_sames[2*J+par](j)=0;
      // 		}
      // 	      ++j;

      // 	    }
      // 	}
      // else if(eigen_values(i).real()>0)
      // 	{
      // 	  Omega_sames[2*J+par](j)=eigen_values(i).real();
      // 	  X_sames[2*J+par].col(j)=eigen_vectors.col(i).head(dim).real();
      // 	  Y_sames[2*J+par].col(j)=eigen_vectors.col(i).tail(dim).real();
      // 	  double norm=X_sames[2*J+par].col(j).squaredNorm()-Y_sames[2*J+par].col(j).squaredNorm();
      // 	  if(norm>0)
      // 	    {
      // 	      X_sames[2*J+par].col(j)/=sqrt(norm);
      // 	      Y_sames[2*J+par].col(j)/=sqrt(norm);
      // 	    }
      // 	  else
      // 	    {
      // 	      Omega_sames[2*J+par](j)=0;
      // 	    }
      // 	  ++j;
      // 	}
      
    }
  
  for(int i=0;i<Lambda.size();i++)
    {
      X_sames[2*J+par].row(i)*=1.0/sqrt(Lambda(i));
      Y_sames[2*J+par].row(i)*=1.0/sqrt(Lambda(i));
    }
  X_sames[2*J+par]=S*X_sames[2*J+par];
  Y_sames[2*J+par]=S*Y_sames[2*J+par];
  return 1;
}


void ERPA::set_AB(int J,int par)
{
  set_ph_np(J,par,ph_nps[2*J+par]);
  set_ph_pn(J,par,ph_pns[2*J+par]);
  int dim_np=ph_nps[2*J+par].size();
  int dim_pn=ph_pns[2*J+par].size();
   
  Anpnps[2*J+par].resize(dim_np,dim_np);
  Apnpns[2*J+par].resize(dim_pn,dim_pn);
  Bnppns[2*J+par].resize(dim_np,dim_pn);

  G_nps[2*J+par].resize(dim_np,dim_np);
  G_pns[2*J+par].resize(dim_pn,dim_pn);

  for(int i=0;i<dim_np;i++)
    {
      for(int j=i;j<dim_np;j++)
	{
	  int bra_p=ph_nps[2*J+par][i].a;
	  int bra_h=ph_nps[2*J+par][i].b;
	  int ket_p=ph_nps[2*J+par][j].a;
	  int ket_h=ph_nps[2*J+par][j].b;
	  G_nps[2*J+par](i,j)=(((bra_p==ket_p)?density_matrix(bra_h,ket_h):0)-((bra_h==ket_h)?density_matrix(bra_p,ket_p):0));
	  G_nps[2*J+par](j,i)=G_nps[2*J+par](i,j);
	}
    }

  for(int i=0;i<dim_pn;i++)
    {
      for(int j=i;j<dim_pn;j++)
	{
	  int bra_p=ph_pns[2*J+par][i].a;
	  int bra_h=ph_pns[2*J+par][i].b;
	  int ket_p=ph_pns[2*J+par][j].a;
	  int ket_h=ph_pns[2*J+par][j].b;
	  G_pns[2*J+par](i,j)=(((bra_p==ket_p)?density_matrix(bra_h,ket_h):0)-((bra_h==ket_h)?density_matrix(bra_p,ket_p):0));
	  G_pns[2*J+par](j,i)=G_pns[2*J+par](i,j);
	}
    }

  
#pragma omp parallel
  {
#pragma omp for
    for(int i=0;i<dim_np;++i)
      {
	for(int j=i;j<dim_np;++j)
	  {
	    int bra_p=ph_nps[2*J+par][i].a;
	    int bra_h=ph_nps[2*J+par][i].b;
	    int ket_p=ph_nps[2*J+par][j].a;
	    int ket_h=ph_nps[2*J+par][j].b;
	    Anpnps[2*J+par](i,j)=a(bra_p,bra_h,ket_p,ket_h,J);
	    Anpnps[2*J+par](j,i)=Anpnps[2*J+par](i,j);
	  }
      }
#pragma omp for
    for(int i=0;i<dim_pn;++i)
      {
	for(int j=i;j<dim_pn;++j)
	  {
	    int bra_p=ph_pns[2*J+par][i].a;
	    int bra_h=ph_pns[2*J+par][i].b;
	    int ket_p=ph_pns[2*J+par][j].a;
	    int ket_h=ph_pns[2*J+par][j].b;
	    Apnpns[2*J+par](i,j)=a(bra_p,bra_h,ket_p,ket_h,J);
	    Apnpns[2*J+par](j,i)=Apnpns[2*J+par](i,j);
	  }
      }
#pragma omp for
    for(int i=0;i<dim_np;++i)
      {
	for(int j=0;j<dim_pn;++j)
	  {
	    int bra_p=ph_nps[2*J+par][i].a;
	    int bra_h=ph_nps[2*J+par][i].b;
	    int ket_p=ph_pns[2*J+par][j].a;
	    int ket_h=ph_pns[2*J+par][j].b;
	    Bnppns[2*J+par](i,j)=b(bra_p,bra_h,ket_p,ket_h,J);
	  }
      }
  }
}

bool ERPA::cal_ce(int J,int par)
{
  set_AB(J,par);
  if(ph_nps[2*J+par].empty()&&ph_pns[2*J+par].empty())
    {
      Anpnps[2*J+par].resize(0,0);
      Apnpns[2*J+par].resize(0,0);
      Bnppns[2*J+par].resize(0,0);
      T_plus_Xs[2*J+par].resize(0,0);
      T_plus_Ys[2*J+par].resize(0,0);
      T_plus_Omegas[2*J+par].resize(0,1);

      T_minus_Xs[2*J+par].resize(0,0);
      T_minus_Ys[2*J+par].resize(0,0);
      T_minus_Omegas[2*J+par].resize(0,1);

      return 0;
    }

  SelfAdjointEigenSolver<MatrixXd> SA_ES1(G_nps[2*J+par]);
  MatrixXd S1=SA_ES1.eigenvectors();
  VectorXd Lambda1=SA_ES1.eigenvalues();

  SelfAdjointEigenSolver<MatrixXd> SA_ES2(G_pns[2*J+par]);
  MatrixXd S2=SA_ES2.eigenvectors();
  VectorXd Lambda2=SA_ES2.eigenvalues();

  Anpnps[2*J+par]=S1.transpose()*Anpnps[2*J+par]*S1;
  Bnppns[2*J+par]=S1.transpose()*Bnppns[2*J+par]*S2;
  for(int i=0;i<Lambda1.size();i++)
    {
      Anpnps[2*J+par].row(i)*=1.0/sqrt(Lambda1(i));
      Anpnps[2*J+par].col(i)*=1.0/sqrt(Lambda1(i));
      Bnppns[2*J+par].row(i)*=1.0/sqrt(Lambda1(i));
    }

  Apnpns[2*J+par]=S2.transpose()*Apnpns[2*J+par]*S2;
  for(int i=0;i<Lambda2.size();i++)
    {
      Apnpns[2*J+par].row(i)*=1.0/sqrt(Lambda2(i));
      Apnpns[2*J+par].col(i)*=1.0/sqrt(Lambda2(i));
      Bnppns[2*J+par].col(i)*=1.0/sqrt(Lambda2(i));
    }

  
  int dim_np=ph_nps[2*J+par].size();
  int dim_pn=ph_pns[2*J+par].size();

  MatrixXcd R(dim_np+dim_pn,dim_np+dim_pn);
  for(int i=0;i<dim_np;++i)
    {
      for(int j=0;j<dim_np;++j)
	{
	  R(i,j)=complex<double>(Anpnps[2*J+par](i,j),0);
	}
      for(int j=0;j<dim_pn;++j)
	{
	  R(i,dim_np+j)=complex<double>(Bnppns[2*J+par](i,j),0);
	}
    }
  for(int i=0;i<dim_pn;++i)
    {
      for(int j=0;j<dim_np;++j)
	{
	  R(dim_np+i,j)=complex<double>(-Bnppns[2*J+par](j,i),0);
	}
      for(int j=0;j<dim_pn;++j)
	{
	  R(dim_np+i,dim_np+j)=complex<double>(-Apnpns[2*J+par](i,j),0);
	}
    }
  ComplexEigenSolver<MatrixXcd> es(R);
  VectorXcd eigen_values=es.eigenvalues();
  MatrixXcd eigen_vectors=es.eigenvectors();
  int T_plus_dim=0,T_minus_dim=0;
  for(int i=0;i<eigen_values.size();i++)
    {
      double temp1=eigen_vectors.col(i).head(dim_np).real().squaredNorm();
      double temp2=eigen_vectors.col(i).tail(dim_pn).real().squaredNorm();
      if(temp1>temp2)
	++T_plus_dim;
      else
	++T_minus_dim;
    }
  T_plus_Omegas[2*J+par].resize(T_plus_dim,1);
  T_minus_Omegas[2*J+par].resize(T_minus_dim,1);

  T_plus_Xs[2*J+par].resize(dim_np,T_plus_dim);
  T_plus_Ys[2*J+par].resize(dim_pn,T_plus_dim);

  T_minus_Ys[2*J+par].resize(dim_np,T_minus_dim);
  T_minus_Xs[2*J+par].resize(dim_pn,T_minus_dim);

  int i_plus=0,i_minus=0;
  for(int i=0;i<eigen_values.size();++i)
    {
      double temp1=eigen_vectors.col(i).head(dim_np).real().squaredNorm();
      double temp2=eigen_vectors.col(i).tail(dim_pn).real().squaredNorm();
      double norm=0;
      if(temp1>temp2)
	{
	  norm=sqrt(temp1-temp2);
	  T_plus_Omegas[2*J+par](i_plus)=eigen_values(i).real();
	  T_plus_Xs[2*J+par].col(i_plus)=eigen_vectors.col(i).head(dim_np).real();
	  T_plus_Ys[2*J+par].col(i_plus)=eigen_vectors.col(i).tail(dim_pn).real();
	  T_plus_Xs[2*J+par].col(i_plus)/=norm;
	  T_plus_Ys[2*J+par].col(i_plus)/=norm;
	  ++i_plus;
	}
      else
	{
	  norm=sqrt(temp2-temp1);
	  T_minus_Omegas[2*J+par](i_minus)=-eigen_values(i).real();
	  T_minus_Ys[2*J+par].col(i_minus)=eigen_vectors.col(i).head(dim_np).real();
	  T_minus_Xs[2*J+par].col(i_minus)=eigen_vectors.col(i).tail(dim_pn).real();
	  T_minus_Ys[2*J+par].col(i_minus)/=norm;
	  T_minus_Xs[2*J+par].col(i_minus)/=norm;
	  ++i_minus;
	}
    }
  for(int i=0;i<Lambda1.size();i++)
    {
      T_plus_Xs[2*J+par].row(i)*=1.0/sqrt(Lambda1(i));
      T_minus_Ys[2*J+par].row(i)*=1.0/sqrt(Lambda1(i));
    }
  for(int i=0;i<Lambda2.size();i++)
    {
      T_plus_Ys[2*J+par].row(i)*=1.0/sqrt(Lambda2(i));
      T_minus_Xs[2*J+par].row(i)*=1.0/sqrt(Lambda2(i));
    }
  T_plus_Xs[2*J+par]=S1*T_plus_Xs[2*J+par];
  T_plus_Ys[2*J+par]=S2*T_plus_Ys[2*J+par];

  T_minus_Ys[2*J+par]=S1*T_minus_Ys[2*J+par];
  T_minus_Xs[2*J+par]=S2*T_minus_Xs[2*J+par];
  return 1;
}
bool ERPA::cal(int type,int J,int par)
{
  return type?cal_ce(J,par):cal_same(J,par);
}

void ERPA::init()
{
  const HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;  
  h_ljjtzs.clear();
  p_ljjtzs.clear();
  for(int i=0;i<=pSystem_Tz_HF->FermiSurface;i++)
    {
      int j=0;
      for(j=0;j<h_ljjtzs.size();j++)
	{
	  if((HFOrbs[i].l==h_ljjtzs[j].l)&&(HFOrbs[i].jj==h_ljjtzs[j].jj)&&(HFOrbs[i].tz==h_ljjtzs[j].tz))
	    {
	      h_ljjtzs[j].ord_list.push_back(i);
	      break;
	    }
	}
      if(j>=h_ljjtzs.size())
	{
	  ljjtz temp;
	  temp.l=HFOrbs[i].l;
	  temp.jj=HFOrbs[i].jj;
	  temp.tz=HFOrbs[i].tz;
	  temp.ord_list.push_back(i);
	  h_ljjtzs.push_back(temp);
	}
    }
  for(int i=pSystem_Tz_HF->FermiSurface+1;i<HFOrbs.size();i++)
    {
      int j=0;
      for(j=0;j<p_ljjtzs.size();j++)
	{
	  if((HFOrbs[i].l==p_ljjtzs[j].l)&&(HFOrbs[i].jj==p_ljjtzs[j].jj)&&(HFOrbs[i].tz==p_ljjtzs[j].tz))
	    {
	      p_ljjtzs[j].ord_list.push_back(i);
	      break;
	    }
	}
      if(j>=p_ljjtzs.size())
	{
	  ljjtz temp;
	  temp.l=HFOrbs[i].l;
	  temp.jj=HFOrbs[i].jj;
	  temp.tz=HFOrbs[i].tz;
	  temp.ord_list.push_back(i);
	  p_ljjtzs.push_back(temp);
	}
    }

  table.resize(HFOrbs.size()); 
  density_matrix_p.resize(p_ljjtzs.size());
  density_matrix_h.resize(h_ljjtzs.size());
  
  for(int i=0;i<h_ljjtzs.size();i++)
    {
      density_matrix_h[i]=MatrixXd::Identity(h_ljjtzs[i].ord_list.size(),h_ljjtzs[i].ord_list.size());
      for(int j=0;j<h_ljjtzs[i].ord_list.size();j++)
	{
	  table[h_ljjtzs[i].ord_list[j]].block_number=i;
	  table[h_ljjtzs[i].ord_list[j]].n=j;
	}
    }
  for(int i=0;i<p_ljjtzs.size();i++)
    {
      density_matrix_p[i]=MatrixXd::Zero(p_ljjtzs[i].ord_list.size(),p_ljjtzs[i].ord_list.size()); 
      for(int j=0;j<p_ljjtzs[i].ord_list.size();j++)
	{
	  table[p_ljjtzs[i].ord_list[j]].block_number=i;
	  table[p_ljjtzs[i].ord_list[j]].n=j;
	}
    }

  
  e_h=MatrixXd::Zero(pSystem_Tz_HF->FermiSurface+1,pSystem_Tz_HF->FermiSurface+1);
  e_p=MatrixXd::Zero(HFOrbs.size()-pSystem_Tz_HF->FermiSurface-1,HFOrbs.size()-pSystem_Tz_HF->FermiSurface-1);
  for(int i=0;i<=pSystem_Tz_HF->FermiSurface;i++)
    {
      e_h(i,i)=HFOrbs[i].e;
    }
  for(int i=pSystem_Tz_HF->FermiSurface+1;i<HFOrbs.size();i++)
    {
      e_p(i-pSystem_Tz_HF->FermiSurface-1,i-pSystem_Tz_HF->FermiSurface-1)=HFOrbs[i].e;
    }
}
void ERPA::cal_density()
{
  for(int i=0;i<density_matrix_p.size();i++)
    {
      density_matrix_p[i]=MatrixXd::Zero(p_ljjtzs[i].ord_list.size(),p_ljjtzs[i].ord_list.size());
    }
  for(int i=0;i<density_matrix_h.size();i++)
    {
      density_matrix_h[i]=MatrixXd::Identity(h_ljjtzs[i].ord_list.size(),h_ljjtzs[i].ord_list.size());
    }

  for(int J=0;J<Jmax;J++)
    for(int par=0;par<2;par++)
      {
	MatrixXd X=G_sames[2*J+par]*X_sames[2*J+par];
	MatrixXd Y=G_sames[2*J+par]*Y_sames[2*J+par];
	MatrixXd X_plus=G_nps[2*J+par]*T_plus_Xs[2*J+par];
	MatrixXd Y_plus=G_pns[2*J+par]*T_plus_Ys[2*J+par];
	MatrixXd X_minus=G_pns[2*J+par]*T_minus_Xs[2*J+par];
	MatrixXd Y_minus=G_nps[2*J+par]*T_minus_Ys[2*J+par];
	
	for(int i=0;i<X.cols();i++)
	  for(int j=0;j<X.cols();j++)
	    {
	      //remove the 1- spurious state
	      if(J==1&&par==1&&(i==0||j==0))
	      	continue;
	      double temp1=0;
	      for(int k=0;k<X.rows();k++)
		{
		  temp1-=X(k,i)*X(k,j);
		}
		  
	      temp1*=0.5;
	      temp1+=((i==j)?1.0:0);
	      
	      for(int k=0;k<X.rows();k++)
		for(int l=0;l<X.rows();l++)
		  {
		    int bra_p=ph_sames[2*J+par][k].a;
		    int ket_p=ph_sames[2*J+par][l].a;

		    int bra_h=ph_sames[2*J+par][k].b;
		    int ket_h=ph_sames[2*J+par][l].b;
		    if((bra_h==ket_h)&&(table[bra_p].block_number==table[ket_p].block_number))
		      {		
			density_matrix_p[table[bra_p].block_number](table[bra_p].n,table[ket_p].n)+=(2*J+1)*temp1*Y(k,i)*Y(l,j)/(p_ljjtzs[table[bra_p].block_number].jj+1);
		      }
		    if((bra_p==ket_p)&&(table[bra_h].block_number==table[ket_h].block_number))
		      {
			density_matrix_h[table[bra_h].block_number](table[bra_h].n,table[ket_h].n)-=(2*J+1)*temp1*Y(k,i)*Y(l,j)/(h_ljjtzs[table[bra_h].block_number].jj+1);
		      }
		  }
	    }


	for(int i=0;i<X_plus.cols();i++)
	  for(int j=0;j<X_plus.cols();j++)
	    {
	      double temp1=0;
	      for(int k=0;k<X_plus.rows();k++)
		{
		  temp1-=X_plus(k,i)*X_plus(k,j);
		}
	      temp1*=0.5;
	      temp1+=((i==j)?1.0:0);

	      for(int k=0;k<Y_plus.rows();k++)
		for(int l=0;l<Y_plus.rows();l++)
		  {
		    int bra_p=ph_pns[2*J+par][k].a;
		    int ket_p=ph_pns[2*J+par][l].a;

		    int bra_h=ph_pns[2*J+par][k].b;
		    int ket_h=ph_pns[2*J+par][l].b;
		    if((bra_h==ket_h)&&(table[bra_p].block_number==table[ket_p].block_number))
		      {
			density_matrix_p[table[bra_p].block_number](table[bra_p].n,table[ket_p].n)+=(2*J+1)*temp1*Y_plus(k,i)*Y_plus(l,j)/(p_ljjtzs[table[bra_p].block_number].jj+1);
		      }
		    if((bra_p==ket_p)&&(table[bra_h].block_number==table[ket_h].block_number))
		      {
			density_matrix_h[table[bra_h].block_number](table[bra_h].n,table[ket_h].n)-=(2*J+1)*temp1*Y_plus(k,i)*Y_plus(l,j)/(h_ljjtzs[table[bra_h].block_number].jj+1);
		      }
		  }
	    }
	    
	for(int i=0;i<X_minus.cols();i++)
	  for(int j=0;j<X_minus.cols();j++)
	    {
	      double temp1=0;
	      for(int k=0;k<X_minus.rows();k++)
		{
		  temp1-=X_minus(k,i)*X_minus(k,j);
		}
	      temp1*=0.5;
	      temp1+=((i==j)?1.0:0);

		 
	      for(int k=0;k<Y_minus.rows();k++)
		for(int l=0;l<Y_minus.rows();l++)
		  {
		    int bra_p=ph_nps[2*J+par][k].a;
		    int ket_p=ph_nps[2*J+par][l].a;

		    int bra_h=ph_nps[2*J+par][k].b;
		    int ket_h=ph_nps[2*J+par][l].b;
		    if((bra_h==ket_h)&&(table[bra_p].block_number==table[ket_p].block_number))
		      {
			density_matrix_p[table[bra_p].block_number](table[bra_p].n,table[ket_p].n)+=(2*J+1)*temp1*Y_minus(k,i)*Y_minus(l,j)/(p_ljjtzs[table[bra_p].block_number].jj+1);
		      }
		    if((bra_p==ket_p)&&(table[bra_h].block_number==table[ket_h].block_number))
		      {
			density_matrix_h[table[bra_h].block_number](table[bra_h].n,table[ket_h].n)-=(2*J+1)*temp1*Y_minus(k,i)*Y_minus(l,j)/(h_ljjtzs[table[bra_h].block_number].jj+1);
		      }
		  }
	    }
      }
}
void ERPA::iteration()
{
  double sigma=0;
  int iter=0;
  do
    {
      vector<MatrixXd> p_old=density_matrix_p;
      vector<MatrixXd> h_old=density_matrix_h;
      #pragma omp parallel for schedule(dynamic) collapse(3)
      for(int J=0;J<Jmax;J++)
	{
	  
	  for(int type=0;type<2;type++)
	    {

	      for(int parity=0;parity<2;parity++)
		{
		  cal(type,J,parity);
		}
	      
	    }
	  
	}
      cal_density();
      set_e();
      sigma=0;
      for(int i=0;i<density_matrix_p.size();i++)
	{
	  sigma+=(density_matrix_p[i]-p_old[i]).norm();
	}
      for(int i=0;i<density_matrix_h.size();i++)
	{
	  sigma+=(density_matrix_h[i]-h_old[i]).norm();
	}
      cout<<"iteration:"<<sigma<<endl;
      iter++;
    }while(sigma>1e-6&&iter<1000);
}

void ERPA::set_q(int J,int Par,int T)
{
  const HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;  
  double Rms=pSystem_Tz_HF->rms();  
  int dim=ph_sames[2*J+Par].size();
  q.resize(dim);
  if(J==1)
    { 
      if(T==0)
	{
	  for(int i=0;i<dim;++i)
	    {
	      q[i]=pSystem_Tz_HF->Q1_s(ph_sames[2*J+Par][i].a,ph_sames[2*J+Par][i].b,Rms);
	    }
	}
      else
	{
	  for(int i=0;i<dim;++i)
	    {
	      q[i]=pSystem_Tz_HF->Q1_v(ph_sames[2*J+Par][i].a,ph_sames[2*J+Par][i].b);
	    }
	}
    }
  else
    {
      for(int i=0;i<dim;++i)
	{
	  q[i]=pSystem_Tz_HF->Q(ph_sames[2*J+Par][i].a,ph_sames[2*J+Par][i].b,J,T);
	}
    }
}
double ERPA::BEJ(int J,int par,int v)
{
  double temp=0;
  int dim=ph_sames[2*J+par].size();
  temp=q.transpose()*G_sames[2*J+par]*(X_sames[2*J+par].col(v)+phase(J)*Y_sames[2*J+par].col(v));
  temp=temp*temp;
  return temp;
}

double ERPA::S(int J,int par,int T)
{
  set_q(J,par,T);
  double temp=0;
  for(int i=0;i<Omega_sames[2*J+par].size();++i)
    {
      temp+=Omega_sames[2*J+par](i)*BEJ(J,par,i);
    }
  return temp;
}

void ERPA::R(int J,int Par,int T,ostream&fout,const double gamma,int points,const double dE)
{
  set_q(J,Par,T);
  double E=0;
  for(int i=0;i<points;++i)
    {
      E=i*dE;
      double temp=0;
      for(int j=0;j<Omega_sames[2*J+Par].size();++j)
	{
	  temp+=0.5*gamma/( (E-Omega_sames[2*J+Par](j))*(E-Omega_sames[2*J+Par](j)) + 0.25*gamma*gamma )/Pi * BEJ(J,Par,j);
	}
      fout<<E<<"\t"<<temp<<endl;
    }
}
